\documentclass{article}
\usepackage{amsmath}
\begin{document}
\section{Shan-Chen Multiple component model}
\label{chapterShanchen} A multi-component system consist of separate
chemical components such as oil and water.Due to their economic
importance, such systems has been studied extensively.\\

Adding a second chemical component without long range interactions
enable us to simulate completely miscible fluid, while completely
immiscible fluid simulation is possible where long range repulsive
interactions are considered in the model.\\

We add a second component by introducing a new index to the relevant
distribution functions: $f^{\sigma}$ represents the distributiong
functions of component $\sigma$. Similar definitions apply to
velocities, densities. To conserve the momentum locally in the
collision step, a composite macroscopic velocity is defined as:

\begin{equation}
u'=\frac{\sum\limits_{\sigma} \frac{1}{\tau_{\sigma}} \sum\limits_a
f_a^{\sigma}
e_a}{\sum\limits_{\sigma}\frac{1}{\tau_{\sigma}}\rho_{\sigma}}
\end{equation}

Uncoupled velocities and densities of individual component $\sigma$
are defined as usual:


\begin{equation}
\rho_a = \sum\limits_{a=0}^8 f_a^{\sigma}
\end{equation}

\begin{equation}
u_{\sigma}=\frac{1}{\rho_{\sigma}} \sum\limits_{a=0}^8 f_a^{\sigma}
e_a
\end{equation}

The force applied on component $\sigma$ is
\begin{equation}
F_{\sigma}(x)=-G \psi_{\sigma}(x,t)\sum\limits_{a}w_a
\psi_{\bar{\sigma}}(x+e_a \Delta t,t)e_a
\end{equation}

Where $\bar{\sigma}$ indicates the other fluid
component.$\psi_{\bar{\sigma}}$ and $\psi_{\sigma}$ are potential
functions and are taken as densities in this example.\\

A repulsive interactive force asks for a positive $G$ value. Large
$G$ value lead to sharp, non-diffusive interface, however give low
numerical stability at the same time.\\

Surface forces are incorporated into the multicomponent model in a
similar way as multi-phase model:

\begin{equation}
F_{surface}^{\sigma} (x,t)= -G_{surface}^{\sigma} \rho(x,t) \sum w_a
S(x+e_a \Delta t,t)e_a
\end{equation}

\subsection{Numerical results}
A series of drops surrounded by another immiscible fluid with same
viscosity are simulated to computed the surface tension. The
simulation results are similar as the SCMP results which confirm the
multi-component model agree with Laplace law.\\

A simulation for two fluid having different viscosities in a channel
is carried out and compared with analytical solutions. An analytical
solution for a sharp interface is :

\begin{equation}
u_x(y)=\begin{cases}
\frac{1}{2\nu_1}G(L^2-a^2)+\frac{1}{2\nu_2}(a^2-y^2)& 0 \le |y| \le a\\
\frac{1}{2\nu_1} G(L^2-y^2) & a \le |y| \le L
\end{cases}
\end{equation}

Where $L$ is the width of the channel and $a$ is the location of
initial interface. The domain is periodic in x direction and bounded
by non-slip walls at $y=0$ and $y=128$ It is initialized with
substance 0 in the middle and substance 1 on both sides near the
walls. Initial density value 1 is applied to substance 0 and 1. The
viscosities of component 0 and 1 are 0.1 and 0.3. Velocity profile
of simulation and analytical solutions is showed in Fig ()
\\%FIGURE!!!!!!!

The overall profile is as expected however it does not give good
agreement with the analytical solution. Density distributions are
plotted in Fig ()%FIGURE!!!!!!!
Concentration decreasing in component 0 and increasing in component
1 are found. This phenomenon is due to the big diffusion on the
interfaces. The diffusivity can be reduced by increasing the
interaction forces however the numerical stability will be lowered
at the same time. Low diffusivity binary simulations are difficult
for Shan-Chen multi-component model due to this numerical stability
problem and deserve future exploration. It's worth mentioning that,
good agreements with analytical solution are obtained for same
viscosity immiscible fluids simulation.


\section{Free Energy Model}
Swift et al (1995) developed a new thermodynamically consistent
binary fluid LB model by introducing an equilibrium state associated
with a free energy functional, corresponding pressure tensors and
chemical potentials. Correct choice of collision rules ensure that
the system evolutes toward to minimization of the free energy
functional.

The thermodynamic properties of a binary fluid system can be
described by a Laudau free energy functional.

\begin{equation}
\Psi=\int_V (\psi_b+\frac{\kappa}{2}(\partial_{\alpha} \phi)^2)dv
+\int_S \psi_s ds
\end{equation}

Where $\psi_b$ is the bulk free energy density and has a form:
\begin{equation}
\psi_b=\frac{c^2}{3} \rho ln\rho +A (-\frac{1}{2}\phi^2+
\frac{1}{4}\phi^4)
\end{equation}

$\phi$ is the order parameter represents the concentration of
components with the definition as:

\begin{equation}
\phi=\frac{\rho_a-\rho_b}{\rho_a+\rho_b}
\end{equation}


This free energy density function can lead to a phase separation
with $\phi=\pm 1$. The gradient term represents an energy
contribution and is related with surface tension. The second
integral is over the surface and describe the interactions between
fluids and the solid surfaces. $\psi_s$ is taken as:
\begin{equation}
\psi_s=-h \phi_s
\end{equation}

$\phi_s$ is the order parameter value at the solid surface.
Minimizing the free energy shows that the gradient in $\phi$ at the
solid surfaces is :

\begin{equation}
\kappa \partial_{\perp} \phi=\frac{d \psi_s}{d \phi_s}=-h
\label{wettingFE}
\end{equation}

The contact angle is related with the parameter $h$ and has a
formulation as:
\begin{equation}
h=\sqrt{2\kappa A} sign(\frac{\pi}{2}-\theta)
\sqrt{cos(\frac{\alpha}{3})[1-cos(\frac{\alpha}{3})]}
\end{equation}

Where $\alpha=cos^{-1}(sin^2 \theta)$.

Pressure tensor and chemical potentials are defined in an usual way:

\begin{eqnarray}
&\mu=\frac{\partial \Psi}{\partial \phi}=A(-\phi+\phi^3)-\kappa
\nabla^2 \phi&\\
&P_{\alpha
\beta}=(P_b-\frac{\kappa}{2}(\partial_{\gamma}\phi)^2-\kappa \phi
\partial_{\gamma \gamma}\phi)\delta_{\alpha \beta}+
\kappa(\partial_{\alpha}\phi)(\partial_{\beta}\phi)& \\
&P_b=\frac{c^2}{3}\rho+A(-\frac{1}{2}\phi^2+\frac{3}{4}\phi^4)&
\end{eqnarray}

The hydrodynamics and thermodynamics of the binary fluids are
described by the Navier-Stokes equations and a convection-diffusion
equation:

\begin{equation}
\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho u)=0
\end{equation}

\begin{equation}
\frac{\partial (\rho u)}{\partial t}+\nabla \cdot (\rho uu)=-\nabla
p_b + \nu \nabla ^2 (\rho u)+ \nabla[\lambda \nabla \cdot (\rho u)]
\end{equation}

\begin{equation}
\frac{\partial \phi}{\partial t}+\nabla (\phi u)= \mu \nabla ^2 \mu
\end{equation}

Where $\lambda= \nu (1-3\frac{c_s^2}{c^2})$.

A new distribution function $g_i(x,t)$ is introduced to describe the
concentration $\phi=\sum\limits_i g_i$ and is related to the
convection and diffusion. The distribution function $f_i(x,t)$ is
related to the fluid density and momentum as usual.

The time evolution equations using MRT scheme can be described as:
\begin{eqnarray}
\text{Collision Step:} & f'_i(x,t)=f_i(x,t)-M^{-1}SM(f_i-f^{eq}_i) & \\
& g'_i(x,t) =g _i^{eq}(x,t)& \\
\text{Streaming Step:} & f_i (x+e_i \Delta t,t+\Delta t) =f'_i(x,t)&
\\
&g_i (x+e_i \Delta t,t+\Delta t) =g'_i(x,t)&
\end{eqnarray}

Appropriate choice of equilibrium distribution fuctions can
reproduce the macroscopic equations in continum limit. The
equilibrium distribution functions should satisfy following
constrains:

\begin{eqnarray}
&\sum\limits_i f_i^{eq} = \rho& \\
&\sum\limits_i f_i^{eq}e_i = \rho u& \\
&\sum\limits_i f_i^{eq}e_i e_i=P + \rho uu & \\
&\sum\limits_i g_i^{eq} = \phi& \\
&\sum\limits_i g_i^{eq}e_i = \phi u& \\
&\sum\limits_i g_i^{eq}e_i e_i=\Gamma \mu I +\phi uu & \\
\end{eqnarray}

\subsection{Wetting boundary condition:}
Bounce Back rule is used to implement the non-slip condition. A
wetting boundary condition is given by equation (\ref{wettingFE})An
appropriate concentration value is assigned to solid nodes to
satisfy the equation(\ref{wettingFE}). $\nabla^2 \phi$ can be
calculated the same way as internal nodes. An example is give in
Figure () %FIGHER!!!!!!!!!!!!!!!!! and algorithm!!!!!!!!!


\subsection{Numerical examples}
A simulation of binary fluids with different viscosity in a channel
is computed with the free energy binary fluids model.

The initial condition is similar as that in Chapter
\ref{chapterShanchen}. The viscosity of component A is 0.01 and 1.2
for component B, the viscosity ratio is 120.

The simulation result (Fig.())%%FIGURE !!!!!!!!!!!!!!!
gives an excellent agreement with the analytical solution. The
concentration profile is also as expected. Two immiscible fluid are
constrains in their respective areas and no unexpected diffusion is
discovered.

We should notice that the interface thickness is around $6$ lattice
units and the interface location is found at $x=$ instead of their
initial position of $x=$ A interface movement of 3 lattice is discovered from the result. %%DATA!!!!!!!!!!!!!!!!!!!!!!!

\section{Color Gradient Model}


An immiscible fluid model developed from Lattice Gas Cellular
Automata was introduced by Gunstensen and Rothman (1991) The
particles in this model are colored either red or blue, therefore it
is normally called color gradient method. The surface tension is
introduced by adding a perturbation to the collision operator while
keep the adherence to the Navier-Stokes equations in homogeous
regions. A recoloring step is involved after surface tension
perturbation calculation in order to achieve zero diffusivity of one
color into the other.

%====================================================================
A general color gradient lattice Boltzmann model can be summarized
as:



We use a D2Q9 lattice model as an example to illustrate the
algorithm of color gradient LBM. An extension to three dimensional
model is straightforward.

We use $f_i^r$, $f_i^b$ and $f_i$ denote the distribution functions
of a red fluid, a blue fluid and their combination respectively. The
densities and velocities can be calculated as:

\begin{eqnarray}
\rho^r=\sum\limits_i f_i^r & \rho^b=\sum\limits_i f_i^b &
\rho=\rho^r+\rho^b \\
&\rho u=\sum\limits_i e_i(f_i^r+f_i^b)&
\end{eqnarray}

An order parameter similar with the free energy model is defined:
\begin{equation}
\phi=\frac{\rho^r-\rho^b}{\rho^r+\rho^b}
\end{equation}

A regular MRT collision step is carried out for $f_i$:
\begin{equation}
f'(x,t)=f(x,t)-M^{-1}SM(f-f^{eq})
\end{equation}

Local kinematic viscosity $\nu$ is computed as a function of the
order parameter $\phi$
\begin{equation}
\nu=\nu_b+\frac{\phi+1}{2}(\nu_r-\nu_b)
\end{equation}

A perturbation is computed to generate the a surface tension.The
surface tension can be expressed as a local anisotropy in the
pressure: the pressure measured normal to the surface is larger than
that tangential to the surface (62) %%CITING!!!!!!!!!!!
The pressure in a LBM is proportional to the density, so surface
tension can be generated by placing preferentially particles in
directions normal to the interface rather than tangential. Mass and
momentum should be conserved.A color gradient G is defined as:

\begin{equation}
G(x,t)=\sum\limits_i e_i \sum\limits_j(f_j^r(x+e_i\Delta
t,t)-f_j^b(x+e_i\Delta t,t))
\end{equation}

Perturbation of the populations are:
\begin{equation}
f''_i(x,t) = f'_i(x,t) +A|G(x,t)|cos(2\theta)
\end{equation}

Where $\theta$ is the angle between the vector $e_i$ and color
gradient $G$. The parameter A is proportional to the magnitude of
surface tension $\sigma$ (D.Rothman and S.Zaleski Lattice Gas
Cellular Automata Cambridge University Press. Cambrige 1997):

\begin{equation}
\sigma=-\frac{192 \rho A}{\lambda}
\end{equation}

A redistribution of color enforce particles to move towards to the
regions occupied by particles with same color. This recoloring step
enable us to achieve the separation of the fluids. It is carried out
by the following maximization problem (Gunstensen and Rothman 1991)

\begin{equation}
W(f_i^{r''},f_i^{b''})=\max\limits_{f_i^{r''},f_i^{b''}}[\sum\limits_i
(f''_{r_i}-f''_{b_i})e_i]
\end{equation}

Subject to constrains:

\begin{eqnarray}
&\rho''_r=\sum\limits_i f_i^{r''}=\rho_r&\\
&f_i^{r''}+f_i^{b''}=f_i&
\end{eqnarray}

This maximization problem makes the color flux $H=\rho_{r}u_r-\rho_b
u_b$ have the same direction of color gradient and conserve the
local mass.

Many algorithms were developed to solve this maximization problem.
We use an algorithm proposed by T\"olke et al (2001). Comparing with
other existing algorithms, it does not generate any velocity
fluctuation for a plain interface and have a good numerical
stability. However the separation of phases is not as strong as the
other schemes.The algorithm is described as:

\begin{equation}
f_i^r=r_i^r+\frac{e_i \cdot G}{|G|} \ast min (f_i^b, f_{i+1}^b,
f_i^r, f_{i+1}^r)
\end{equation}

\begin{equation}
f_{i+1}^r=r_{i+1}^r-\frac{e_i \cdot G}{|G|} \ast min (f_i^b,
f_{i+1}^b, f_i^r, f_{i+1}^r)
\end{equation}

\begin{equation}
f_i^b=r_i^b-\frac{e_i \cdot G}{|G|} \ast min (f_i^b, f_{i+1}^b,
f_i^r, f_{i+1}^r)
\end{equation}

\begin{equation}
f_{i+1}^b=r_{i+1}^b+\frac{e_i \cdot G}{|G|} \ast min (f_i^b,
f_{i+1}^b, f_i^r, f_{i+1}^r)
\end{equation}

A general color gradient lattice Boltzmann model can be summarized
as:

\begin{itemize}
\item Single phase collision.
    %\begin{equation}
    %f_i=f_i^r+f_i^b
    %\end{equation}
\item Add a surface tension perturbation to $f'_i$ obtaining $f''_i$
\item Recoloring
\item Steaming
\end{itemize}

\section{Verification}
\subsection{Poiseuille flow-Color Gradient Model}

Good agreement with analytical solution is obtained as showed in
Fig()%%FIGURE!!!!!!!!!!!!!!!!
The maximum viscosity ratio for color gradient model is around 120
(Static interface) which is similar to the free energy model.
However, the interface shifting observed in color gradient model is
around 1.5 lattice units which is significantly smaller than the
free energy model. Concentration profile illustrate the diffusive
properties of color gradient (Fig ) %FIGURE!!!!!!!!!!!!!
The recoloring algorithm separates two immiscible fluids very well
and gives 0 diffusivity as we expected.

\subsection{Capillary filling dynamics}
Since Shan-Chen multi-component model cannot deal with binary fluids
with different viscosities, the following comparisons are carried
out between color gradient model and free energy model.

A capillary tube filled with binary fluids with viscosity contrast
is simulated. A hydrophilic inner surface is associated with the
tube. The energy liberated in wetting the surface is used to drive
the fluid up the tube.

An analysis of capillary penetration is given by Pooley et al
(2009). The length of the liquid column that penetrates the
capillary could be calculated by surface tension $\sigma_{lg}$,
capillary tube width $H$, dynamic contact angle $\theta$ and liquid
viscosity $\eta$:

\begin{equation}
l=(\frac{\sigma_{lg}H cos
\theta}{3\eta})^{\frac{1}{2}}(t+t_0)^{\frac{1}{2}}
\label{Washburnequation}
\end{equation}

A system of $512 \times 32$ lattice units with periodic boundary
condition in x direction is used in verification. Non-slip wetting
boundary conditions are involved at the upper and lower sides of the
tube.The simulation setup is showed in Fig( ) %FIGURE!!!!!!!!!!!
An equilibrium gradient at the boundary $\partial_n
\rho|_{bc}=-0.144$ is chosen for free energy model such that the
equilibrium contact angle is $\theta=60^\circ$. The interface is
initialized at $l=50$. Capillary filling distance are compared
against with the Washburn equation (\ref{Washburnequation}).

%%FIGURE!!!!!!!!!!!!!

The distance of binary fluids along with capillary as a function of
time. circle: $\nu_l=$, $\nu_g=$. square solid and dashed line are
theoretical predictions.

We should notice that the maximum viscosity ratio of both free
energy model and color gradient model decrease from around 120
(static interface) to around 20 for dynamical interface simulation.

\subsection{Fingering}
Fingering phenomenons was studied since 1960s due to its relation with the recovery of oil and the transport of substance in porous media.We investigate it by lattice Boltzmann binary models. One
fluids is pushed by another one, with different viscosities, along a
channel with non-slip wall. A body force $G$ is applied on the
fluids. A growing finger of the driving fluid. If the capillary
number $Ca$ is big enough.

\begin{equation}
Ca=\frac{u_t \rho \nu_2}{\sigma}
\end{equation}

Where $u_t$ is the velocity of the tip of the finger, $\nu_2$ is the
viscosity of driving fluid. $\sigma$ is the surface tension.

\end{document}
